fs = 100000;    % 采样频率为 100 kHz
wp = [5000 15000]*2/fs; % 通带归一化截止频率
ws = [4000 16000]*2/fs; % 阻带归一化截止频率
[n,Wn] = buttord(wp, ws, 3, 60); % 计算滤波器阶数和归一化截止频率
[b,a] = butter(n, Wn, 'bandpass'); % 计算滤波器系数

f = logspace(1, 5); % 生成从10^1到10^5的对数均匀分布的100个数字，单位为赫兹
w = f*2*pi; % 将频率转换为角频率

h = freqs(b, a, w); % 求取滤波器在频域的复频率响应
mag = abs(h); % 复频率响应的幅度
phase = angle(h); % 复频率响应的相位
phasedeg = phase*180/pi; % 相位转换为角度

subplot(2,1,1); % 将画图区域分成两个子图，第一个子图
semilogx(f, 20*log10(mag)); % 使用对数坐标画幅频曲线
grid on; % 添加网格线
xlabel('Frequency (Hz)'); % 设置x轴标签
ylabel('Magnitude (dB)'); % 设置y轴标签

subplot(2,1,2); % 第二个子图
semilogx(f, phasedeg); % 使用对数坐标画相频曲线
grid on; % 添加网格线
xlabel('Frequency (Hz)'); % 设置x轴标签
ylabel('Phase (degrees)'); % 设置y轴标签
